Woodland Restoration Survey

Load packages & data

library(tidyverse)
library(sf)
library(vegan)
library(goeveg) # needed for dimcheckMDS to create scree plot.
library(here)
library(pairwiseAdonis)
library(RVAideMemoire)

`%!in%` <- Negate(`%in%`)

source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
# source(here("Scripts", "4_Covariate_Data_Prep.R"))
# source(here("Scripts", "5_Prep_Species.R"))

Load WRS data

wrs <- read.csv(here("Data", "Wildlife", "Raw_bird_survey_data", "2022", "WRS_Data_20220816.csv"))

Clean data

wrs <- wrs %>%
  filter(Temp_issue != 1 | is.na(Temp_issue)) %>% # remove rows that need to be fixed
  filter(Species != "") %>% # remove blank species (temporary)
  filter(Flyover != 1 | is.na(Flyover)) %>% # remove flyovers
  filter(!Initial_detection_time_interval %in% c("After count", "Before count")) %>% # remove detection occurring before/after count
  filter(!Initial_detection_distance_bin %in% c("Beyond woodlot"))

# Impute missing distance bins from best guess
wrs <- wrs %>%
  mutate(Initial_detection_distance_bin = if_else(Initial_detection_distance_bin == "", ">50 m but in woodlot", Initial_detection_distance_bin))

# Split species code and common name
wrs <- wrs %>%
  separate(Species, c("Sp", "Common_name"), sep = ": ")

# Drop unidentified birds
wrs <- wrs[!grepl("Unidentified", wrs$Common_name), ]

# Get unique sites (create standalone df)
wrs_sites <- unique(wrs$Site_ID_full) %>%
  as.data.frame() %>%
  rename("Site_ID_full" = ".")

# Get unique surveys (add to main df and create standalone)
wrs <- wrs %>%
  mutate(Survey_ID_full = paste0(Site_ID_full, " - ", Date))

wrs_surveys <- unique(wrs$Survey_ID_full) %>%
  as.data.frame() %>%
  rename("Survey_ID_full" = ".") %>%
  separate(col = Survey_ID_full, remove = FALSE, sep = " - ", into = c("SurveyType", "Park", "Site", "Date")) %>%
  mutate(Site_ID_full = paste(SurveyType, Park, Site, sep = " - ")) %>%
  dplyr::select(Survey_ID_full, Site_ID_full)

# TODO add survey covariates to survey table (wind, time, etc)

# Join site info from shapefile
# Make site IDs in site info df match those in bird data
pts <- st_read(here("Projects", "Woodland_Bird_Surveys", "Survey_locations_attributed.shp")) %>%
  mutate(Park = ifelse(Park == "Cleary Lake", "Cleary", Park))
site_info <- pts

# Drop the truncated park names from part after final "-"
site_info$Site_ID_full <- gsub("CH", "", site_info$SurveyID)
site_info$Site_ID_full <- gsub("BK", "", site_info$Site_ID_full)
site_info$Site_ID_full <- gsub("LR", "", site_info$Site_ID_full)

# Remove non-trailing zeros
site_info$Site_ID_full <- gsub("(?<![0-9])0+", "", site_info$Site_ID_full, perl = TRUE)

# Tweak survey type names in full ID to match
site_info$Site_ID_full <- gsub("WRS Point Count", "Woodland Restoration Survey", site_info$Site_ID_full)
site_info$Site_ID_full <- gsub("Point Count MB", "Mountain Bike Trail Survey", site_info$Site_ID_full)

site_info <- site_info %>%
  dplyr::select(Site_ID_full, Type, YearPlant)

# Manually append info on points added later
# TODO add these to actual shapefile and confirm values
new_site_data <- data.frame(
  Site_ID_full = c("Woodland Point Count - Baker - 1", "Woodland Point Count - Baker - 2"),
  Type = c("Mature forest", "Mature forest"),
  YearPlant = c(1953, 1953)
)
site_info <- bind_rows(site_info, new_site_data)

wrs_sites <- left_join(wrs_sites, site_info, by = join_by(Site_ID_full))

# Arrange data
wrs_surveys <- left_join(wrs_surveys, site_info, by = join_by(Site_ID_full)) %>%
  arrange(Survey_ID_full) %>%
  dplyr::select(-c(Survey_ID_full, Site_ID_full, geometry)) %>%
  mutate(
    Time = 2022 - YearPlant,
    Type = as.factor(Type)
  ) %>%
  dplyr::select(Time, Type)

wrs_surveys$Type <- factor(wrs_surveys$Type, levels = c("Mature forest", "LBR", "Seedling", "Natural succession")) # reorder factor levels, moving mature first

wrs_species <- wrs %>%
  group_by(Survey_ID_full, Sp) %>%
  summarize(Count = n()) %>%
  spread(Sp, Count, fill = 0) %>%
  arrange(Survey_ID_full) %>%
  ungroup() %>%
  dplyr::select(-Survey_ID_full)

Fit NDMS

# Generate distance matrix with default settings (bray-curtis)
dist_mat <- vegdist(wrs_species)

# Determine best dimensionailty by comparing stress levels with scree plot
screeplot <- dimcheckMDS(wrs_species,
  distance = "bray", # distance measure used for metaMDS analysis
  trymax = 200,
  k = 6
) # maximum number of dimensions, default is 6. could in theory test up to max number of species.

Stress with just 2 dimensions is >0.2, which is unusable according to some. Going to use 3 dimensions.

set.seed(1)
nmds3 <- metaMDS(
  comm = wrs_species, # Define the data to use
  distance = "bray", # Specify distance
  k = 3,
  try = 1000, trymax = 1000, # Number of random starts
  autotransform = TRUE, # allow default for community data
  noshare = TRUE, # allow default for community data
  wascores = TRUE
) # allow default for community data

nmds3

Evaluate NDMS

nmds3$ndim # how many nMDS axes were created?

## [1] 3

nmds3$converged # did it find a convergent solution?

## [1] 11

# nmds3$points # scores for each location (as samples)
# nmds3$species # scores for each animal order (as variables)

A Shepard diagram is a scatterplot of distances between points in an nMDS final configuration (observed dissimilarity) against the original dissimilarities (ordination distance). For a nMDS, the line of best fit produced is monotonic, in that it can only step up or remain constant. The points in a good Shepard diagram should be a clean curve or straight line. Large scatter around the line suggests that original dissimilarities are not well preserved in the reduced number of dimensions. If your plot resembles a step-wise or inverse L-shaped function, it would be wise to not interpret this data and instead try to solve for a degenerate solution (not covered in this course). However, please note if you have a small data set (not many observations) your data may look step-wise in a Shepard diagram. Note, you will need to produce a new Shepard diagram when altering the number of dimensions being plotted. Confusingly, the Shepard diagram is called stressplot in R!

stressplot(nmds3)

TODO: Looks like I might have a step. How problematic is this and what to do about it?

What happens if I run on version without species only observed once?

sp_singletons <- wrs %>%
  group_by(Sp) %>%
  summarize(n = n()) %>%
  filter(n == 1) %>%
  pull(Sp)

wrs_species_no_singletons <- wrs %>%
  filter(Sp %!in% sp_singletons) %>%
  group_by(Survey_ID_full, Sp) %>%
  summarize(Count = n()) %>%
  spread(Sp, Count, fill = 0) %>%
  arrange(Survey_ID_full) %>%
  ungroup() %>%
  dplyr::select(-Survey_ID_full)

screeplot_no_singletons <- dimcheckMDS(wrs_species_no_singletons,
  distance = "bray", # distance measure used for metaMDS analysis
  trymax = 200,
  k = 6
) # maximum number of dimensions, default is 6. could in theory test up to max number of species.

# screeplot_no_singletons

# Stress from just 2 dimensions is still >0.2. going to go with 3 dimensions.
nmds_no_singletons <- metaMDS(
  comm = wrs_species_no_singletons, # Define the data to use
  distance = "bray", # Specify distance
  k = 3,
  try = 1000, trymax = 1000, # Number of random starts
  autotransform = TRUE, # allow default for community data
  noshare = TRUE, # allow default for community data
  wascores = TRUE
) # allow default for community data

# nmds_no_singletons

stressplot(nmds_no_singletons)

Eliminating singletons didn’t seem to change the Shepard diagram much, so I’m going to use full dataset for rest of analysis.

Plotting ordination

Basic biplot with ellipses around all sites of each type.

# Define colors for treatments
colvec <- c("green4", "orange", "maroon", "royalblue")
scl <- 3 # scaling term

# Generate biplot with ellipses
plot(nmds3, type = "n", las = 1)
text(nmds3, display = "species", cex = 0.55, col = "grey40")
points(nmds3,
  display = "sites", col = colvec,
  scaling = scl, pch = 21, cex = 0.75, bg = colvec
)
ordiellipse(nmds3, wrs_surveys$Type,
  display = "sites",
  kind = "ehull", label = F,
  lwd = 2, col = colvec
)
legend("topright",
  legend = levels(wrs_surveys$Type), bty = "n",
  col = colvec, pch = 21, pt.bg = colvec, cex = .5
)

Biplot where ellipses represent error around centroid: “When kind=‘sd’ or ‘se’, the ellipse represents a calculated region of error around each group’s centroid. The confidence level is set with conf= and can be understood as how single-dimensional error bars or confidence intervals might be used in univariate plots to estimate the precision of estimated mean or centroid (se) or connect variation in the sampled data to the spread of the broader population (sd)” (Source)

plot(nmds3, type = "n", las = 1)
ordiellipse(nmds3, wrs_surveys$Type,
  display = "sites",
  kind = "se", conf = 0.95,
  label = F, col = colvec,
  draw = "polygon", alpha = 50, scaling = scl
)
points(nmds3,
  display = "sites", col = colvec,
  scaling = scl, pch = 21, cex = 0.75, bg = colvec
)
text(nmds3, display = "species", scaling = scl, cex = 0.55, col = "grey40")
legend("topright",
  legend = levels(wrs_surveys$Type), bty = "n",
  col = colvec, pch = 21, pt.bg = colvec, cex = .5
)

Now generate same plot but with filtered species for clarity:

plot(nmds_no_singletons, type = "n", las = 1)
ordiellipse(nmds_no_singletons, wrs_surveys$Type,
  display = "sites",
  kind = "se", conf = 0.95,
  label = F, col = colvec,
  draw = "polygon", alpha = 50, scaling = scl
)
points(nmds_no_singletons,
  display = "sites", col = colvec,
  scaling = scl, pch = 21, cex = 0.75, bg = colvec
)
text(nmds_no_singletons, display = "species", scaling = scl, cex = 0.55, col = "grey40")
legend("topright",
  legend = levels(wrs_surveys$Type), bty = "n",
  col = colvec, pch = 21, pt.bg = colvec, cex = .5
)

Mature forest separates from restored sites (aligned with species like WOTH, EAWP, HAWO, SCTA, OVEN). Restoration types have overlap here, so not major differences (but doesn’t account for site age, etc.).

Post hoc tests

“The key function for testing groups (and gradients) against an ordination ad hoc–after the ordination has been fit and the latent variable described–is envfit. It is set up similarly to the ordi* functions, but since it is a test, it uses formula notation.”

envfit(nmds3 ~ wrs_surveys$Type + wrs_surveys$Time, choices = c(1:3)) # use all 3 dimensions when hypothesis testing (default is just 2)

## 
## ***VECTORS
## 
##                      NMDS1     NMDS2     NMDS3     r2 Pr(>r)    
## wrs_surveys$Time -0.851510 -0.038286  0.522940 0.4408  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                                      NMDS1   NMDS2   NMDS3
## wrs_surveys$TypeMature forest      -0.4078 -0.0574  0.1531
## wrs_surveys$TypeLBR                 0.0404 -0.1170 -0.2384
## wrs_surveys$TypeSeedling            0.2338 -0.0248 -0.0442
## wrs_surveys$TypeNatural succession  0.1155  0.1502  0.0077
## 
## Goodness of fit:
##                      r2 Pr(>r)    
## wrs_surveys$Type 0.2257  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Test indicates that the avian community is significantly different in at least 1 of the woodland types. Time since revegetation also important. Note that envfit tests each factor independently.

We now check pairwise comparisons to see which woodland types differ.

pairwise.factorfit(nmds3, wrs_surveys$Type)

## 
##  Pairwise comparisons using factor fitting to an ordination 
## 
## data:  nmds3 by wrs_surveys$Type
## 999 permutations 
## 
##                    Mature forest LBR   Seedling
## LBR                0.004         -     -       
## Seedling           0.003         0.170 -       
## Natural succession 0.003         0.092 0.092   
## 
## P value adjustment method: fdr

Mature forest differs significantly from all other types. LBR and seedling not significantly different. Some suggestion that both LBR and Seedling sites differ from Natural Succession sites, but not significant at p < 0.05. (p = 0.084).

It would be good to analyze interaction between time and type. I think adonis can be used for this.

adonis2(wrs_species ~ wrs_surveys$Type * wrs_surveys$Time, by = "terms")

Type is highly significant (communities in at least one forest type are different from others). Time alone isn’t significant after controlling for type, but there is an interaction effect (influence of time is mediated by type). Suggests different trajectories for different types…

I don’t fully understand this command, but by = “onedf” has interesting output. From help files “by =”terms” will assess significance for each term (sequentially from first to last), setting by = “margin” will assess the marginal effects of the terms (each marginal term analysed in a model with all other variables), by = “onedf” will analyse one-degree-of-freedom contrasts sequentially, by = NULL will assess the overall significance of all terms together. The arguments is passed on to anova.cca.”

adonis2(wrs_species ~ Type * Time, data = wrs_surveys, by = "onedf")

This is kinda interesting. Suggests significant interaction with time and seedling but not LBR (because LBR jump starts process and less variability over time?). Also LBR not signifcantly diff.

Pairwise comparisons:

pairwise.adonis2(wrs_species ~ Type, data = wrs_surveys)

## $parent_call
## [1] "wrs_species ~ Type , strata = Null , permutations 999"
## 
## $`Natural succession_vs_Mature forest`
##          Df SumOfSqs      R2      F Pr(>F)    
## Type      1   1.7598 0.11425 6.9651  0.001 ***
## Residual 54  13.6440 0.88575                  
## Total    55  15.4039 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Natural succession_vs_Seedling`
##          Df SumOfSqs      R2      F Pr(>F)  
## Type      1   0.5021 0.03035 1.8778  0.028 *
## Residual 60  16.0431 0.96965                
## Total    61  16.5452 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Natural succession_vs_LBR`
##          Df SumOfSqs     R2      F Pr(>F)  
## Type      1    0.431 0.0412 1.6329  0.041 *
## Residual 38   10.029 0.9588                
## Total    39   10.460 1.0000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Mature forest_vs_Seedling`
##          Df SumOfSqs      R2      F Pr(>F)    
## Type      1   2.4625 0.13554 9.7208  0.001 ***
## Residual 62  15.7057 0.86446                  
## Total    63  18.1681 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Mature forest_vs_LBR`
##          Df SumOfSqs      R2      F Pr(>F)    
## Type      1   1.1478 0.10589 4.7372  0.001 ***
## Residual 40   9.6919 0.89411                  
## Total    41  10.8397 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Seedling_vs_LBR
##          Df SumOfSqs      R2      F Pr(>F)
## Type      1   0.3466 0.02787 1.3187  0.175
## Residual 46  12.0909 0.97213              
## Total    47  12.4376 1.00000              
## 
## attr(,"class")
## [1] "pwadstrata" "list"

pairwise.adonis2(wrs_species ~ Type * Time, data = wrs_surveys)

## $parent_call
## [1] "wrs_species ~ Type * Time , strata = Null , permutations 999"
## 
## $`Natural succession_vs_Mature forest`
##          Df SumOfSqs      R2      F Pr(>F)    
## Type      1   1.7598 0.11425 6.9163  0.001 ***
## Time      1   0.1583 0.01028 0.6221  0.844    
## Residual 53  13.4857 0.87548                  
## Total    55  15.4039 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Natural succession_vs_Seedling`
##           Df SumOfSqs      R2      F Pr(>F)  
## Type       1   0.5021 0.03035 1.9138  0.027 *
## Time       1   0.2610 0.01578 0.9949  0.473  
## Type:Time  1   0.5655 0.03418 2.1555  0.011 *
## Residual  58  15.2166 0.91970                
## Total     61  16.5452 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Natural succession_vs_LBR`
##           Df SumOfSqs      R2      F Pr(>F)  
## Type       1   0.4310 0.04120 1.6146  0.056 .
## Time       1   0.1529 0.01462 0.5730  0.902  
## Type:Time  1   0.2678 0.02560 1.0032  0.427  
## Residual  36   9.6085 0.91858                
## Total     39  10.4602 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Mature forest_vs_Seedling`
##          Df SumOfSqs      R2      F Pr(>F)    
## Type      1   2.4625 0.13554 9.9890  0.001 ***
## Time      1   0.6682 0.03678 2.7107  0.003 ** 
## Residual 61  15.0375 0.82768                  
## Total    63  18.1681 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Mature forest_vs_LBR`
##          Df SumOfSqs      R2      F Pr(>F)    
## Type      1   1.1478 0.10589 4.7473  0.001 ***
## Time      1   0.2624 0.02421 1.0853  0.367    
## Residual 39   9.4295 0.86990                  
## Total    41  10.8397 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Seedling_vs_LBR
##           Df SumOfSqs      R2      F Pr(>F)  
## Type       1   0.3466 0.02787 1.3666  0.148  
## Time       1   0.5011 0.04029 1.9757  0.021 *
## Type:Time  1   0.4295 0.03453 1.6933  0.045 *
## Residual  44  11.1603 0.89731                
## Total     47  12.4376 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Restoration sites vs reference: Mature forest differs from all other types. Bird communities in restoration sites have not achieved a composition similar to the mature reference sites.

Active vs. passive restoration: Natural succession and seedling sites significantly different. Difference between LBR and natural succession sites is marginally significant. There is a significant interaction between type and time when comparing natural succession and seedling sites (community changes at different rates in seedling sites than natural succession sites), but not for LBR.

Large bareroot vs seedling plantings: No significant difference based on restoration type alone, but the significant interaction between type and time suggests restoration type does effect the rate at which species composition changes.

Check for dispersion

From adnois help file: “Anderson (2001, Fig. 4) warns that the method may confound location and dispersion effects: significant differences may be caused by different within-group variation (dispersion) instead of different mean values of the groups (see Warton et al. 2012 for a general analysis). However, it seems that adonis2 is less sensitive to dispersion effects than some of its alternatives (anosim, mrpp). Function betadisper is a sister function to adonis2 to study the differences in dispersion within the same geometric framework.”

Also see end of this help vid: https://www.youtube.com/watch?v=oLf0EpMJ4yA

bd <- betadisper(vegdist(wrs_species), wrs_surveys$Type)
anova(bd)

permutest(bd)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups      3 0.02460 0.0081993 0.9149    999  0.416
## Residuals 100 0.89621 0.0089621

Doesn’t seem to be different within-group variation, which is good, I think.

Visulize distances

Each point in these plots is a paired comparison between each site and each mature forest sites. Site pairs with lower Bray-Curtis differences have more similar bird communities (less difference). A values of 1 indicates that the sites have no species in common; a value of 0 indcates that the sites have the same species composition.

# Jitter plot. All woodland types vs mature forest pots.
dist_mat_tidy %>%
  filter(str_detect(comparison, "Mature")) %>%
  ggplot(aes(x = comparison, y = value)) +
  geom_jitter(aes(color = as.numeric(time_difference)), size = .5, width = 0.4, alpha = .9) +
  stat_summary(
    fun.data = median_hilow, color = "red", size = 1,
    fun.args = list(conf.int = 0.50)
  ) +
  labs(x = NULL, y = "Bray-Curtis differences", color = "Time difference between sites (years)") +
  # scale_x_discrete(breaks=c("early", "late", "same"),
  #                  labels=c("Inter-mouse\ndistances at\n9 dpw",
  #                           "Inter-mouse\ndistances at\n141 dpw",
  #                           "Intra-mouse\ndistances between\n9 and 141 dpw")) +
  scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
  theme_classic() +
  scale_color_viridis_c() +
  theme(axis.text.x = element_text(angle = 90))

The mature forest vs. mature forest comparison is a measurement of normal community variability within reference sites (even mature forest sites differ form one another to a certain degree). On average the LBR sites are more similar to mature forests than the seedling and natural succession sites, but the effect is small and significance not tested.

Similar plot but instead of coding time differences with color we plot time difference on the x axis. Time difference is the difference in site ages in years, ie how much older the mature forests are than the restoration).

# Time difference on x-axis. Tall facets.
dist_mat_tidy %>%
  filter(str_detect(comparison, "Mature")) %>%
  ggplot(aes(x = time_difference, y = value)) +
  geom_point(aes(color = Type.y), size = .9, alpha = .25) +
  geom_smooth(method = lm, se = TRUE, color = "black") +
  stat_summary(
    fun.data = median_hilow, color = "grey30", size = .3,
    fun.args = list(conf.int = 0.50)
  ) +
  # stat_summary(fun.data=median_hilow, color="red", size=1, fun.args = list(conf.int=0.50)) +
  labs(x = NULL, y = "Bray-Curtis differences") +
  scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
  scale_x_reverse() +
  facet_wrap(~Type.y, nrow = 1) +
  ggthemes::theme_clean() +
  labs(
    x = "Time difference (mature forest head-start)",
    y = "Community dissimilarity (Bray-Curtis distance)",
    color = "Woodland plot type\ncompared to\nmature forest plot"
  )

# Time difference on x-axis. Long facets.
dist_mat_tidy %>%
  filter(str_detect(comparison, "Mature")) %>%
  ggplot(aes(x = time_difference, y = value)) +
  geom_point(aes(color = Type.y), size = .9, alpha = .25) +
  geom_smooth(method = lm, se = TRUE, color = "black") +
  stat_summary(
    fun.data = median_hilow, color = "grey30", size = .3,
    fun.args = list(conf.int = 0.50)
  ) +
  scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
  scale_x_reverse() +
  facet_wrap(~Type.y, nrow = 4) +
  ggthemes::theme_clean() +
  labs(
    x = "Time difference (mature forest head-start)",
    y = "Plot community dissimilarity (vs. mature forest plots)",
    color = "Woodland plot type\ncompared to\nmature forest plot"
  )

# Time on x-axis. Overlayed. With "cohorts".
dist_mat_tidy %>%
  filter(str_detect(comparison, "Mature")) %>%
  ggplot(aes(x = time_difference, y = value)) +
  geom_point(aes(color = Type.y), size = .9, alpha = .1) +
  geom_smooth(aes(group = Type.y, color = Type.y), method = lm, se = FALSE, size = 2) +
  stat_summary(
    fun.data = median_hilow, aes(fill = Type.y), color = "grey20", size = .4, pch = 24, stroke = .5,
    fun.args = list(conf.int = 0.50)
  ) +
  labs(x = NULL, y = "Bray-Curtis differences") +
  scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
  scale_x_reverse() +
  theme_classic() +
  labs(
    x = "Time difference (mature forest head-start)",
    y = "Plot community dissimilarity (vs. mature forest plots)"
  )

# Time on x-axis. Overlayed. Without "cohorts".
dist_mat_tidy %>%
  filter(str_detect(comparison, "Mature")) %>%
  ggplot(aes(x = time_difference, y = value)) +
  geom_point(aes(color = Type.y), size = 1.75, alpha = .1) +
  geom_smooth(aes(group = Type.y, color = Type.y), method = lm, se = TRUE, size = 1.25) +
  geom_smooth(method = lm, se = FALSE, color = "grey40", size = .5) +
  labs(x = NULL, y = "Bray-Curtis differences") +
  scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
  scale_x_reverse() +
  theme_classic() +
  labs(
    x = "Time difference (mature forest head-start)",
    y = "Plot community dissimilarity (vs. mature forest plots)",
    color = "Woodland plot type\ncompared to\nmature forest plot"
  )

These plots show that restoration sites and natural succession sites are become more similar to mature forest as they age. LBR sites start out more similar to mature forests than seedling sites. It takes about 20 years for the seedling sites to catch up with LBR sites. Natural succession community progression towards mature sites is slower than LBR and seedling.